Feature selection analysis

Reproducibility package for the paper …

R. Torkar and R. Berntsson Svensson

2020-05-28

Introduction

First prepare the data, check for NAs and look at some descriptive statistics. Since we’re using Excel we should be very careful when loading the data to see that nothing goes wrong (data manipulated in an arbitrary way, uncompatible data types, etc.)

d <- read.xlsx("Features.xlsx", sheet = "Features")

This is how the data looks like. In short, \(11110\) rows and \(10\) variables.

str(d)
## 'data.frame':    11110 obs. of  10 variables:
##  $ ID                    : num  1 2 3 4 5 6 7 8 9 10 ...
##  $ State                 : chr  "Elicited, Prio, Dropped" "Elicited, Prio, Dropped" "Elicited, Prio, Planned, Implemented, Tested, Released" "Elicited, Prio, Dropped" ...
##  $ Team.priority         : num  88 83 957 79 88 76 74 73 72 71 ...
##  $ Critical.feature      : chr  "No" "No" "No" "No" ...
##  $ Business.value        : chr  "No value" "No value" "No value" "No value" ...
##  $ Customer.value        : chr  "No value" "No value" "No value" "No value" ...
##  $ Stakeholders          : num  1 1 0 0 0 1 1 1 1 1 ...
##  $ Key.customers         : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Dependency            : chr  "No" "No" "No" "No" ...
##  $ Architects.involvement: chr  "None" "None" "None" "None" ...

Each row has a unique ID. The State, which is our outcome (dependent variable), shows how far the feature got in the process. It is of an ordered categorical type, which should be modeled with a cumulative or adjacent category likelihood. There are seven categories:

  1. Elicited, Dropped
  2. Elicited, Prio, Dropped
  3. Elicited, Prio, Planned, Dropped
  4. Elicited, Prio, Planned, Implemented, Dropped
  5. Elicited, Prio, Planned, Implemented, Tested, Dropped
  6. Elicited, Prio, Planned, Implemented, Tested, Released

Team.priority is the relative priority the feature got, \(\mathbb{N} = \{0,\ldots,1000\}\). Critical.feature is a simple ‘Yes’/‘No’ answer (\(\mathbb{Z}_2\)). Business.value and Customer.value are also ordered categorical with three levels and a fourth level called ‘No value’:

  1. No value
  2. Valuable
  3. Important
  4. Critical

Stakeholders have integers, i.e., \(\mathbb{N} = \{0,\ldots,10\}\), and Key.customers the same, but with a different set, i.e., \(\mathbb{N} = \{0,\ldots,60\}\).

Finally, Dependency is \(\mathbb{Z}_2\), while Architects.involvement is ordered categorical:

  1. None
  2. Simple
  3. Monitoring
  4. Active Participation
  5. Joint Design

All ordered categorical predictors (independent variables) can be modeled as monotonic or category-specific effects, if necessary [1[1] P.-C. Bürkner and E. Charpentier, “Modelling monotonic effects of ordinal predictors in Bayesian regression models,” British Journal of Mathematical and Statistical Psychology, doi: 10.1111/bmsp.12195.].

table(is.na(d))
## 
##  FALSE 
## 111100

No NAs in the dataset. However, that doesn’t mean that we don’t have NAs. Some of the coding can be a representation of NA, e.g., ‘No value’. In this particular case we know that ‘No value’ and ‘None’ in the dataset actually are values and not a representation of NA.

Finally, we should set correct data types on all predictors.

# ordered categorical
d$State <- factor(d$State, 
                  levels = c("Elicited, Dropped", 
                             "Elicited, Prio, Dropped", 
                             "Elicited, Prio, Planned, Dropped", 
                             "Elicited, Prio, Planned, Implemented, Dropped", 
                             "Elicited, Prio, Planned, Implemented, Tested, Dropped", 
                             "Elicited, Prio, Planned, Implemented, Tested, Released"), 
                  ordered = TRUE)

d$Business.value <- factor(d$Business.value, 
                           levels = c("No value",
                                      "Valuable",
                                      "Important",
                                      "Critical"), 
                           ordered = TRUE)

d$Customer.value <- factor(d$Customer.value, 
                           levels = c("No value",
                                      "Valuable",
                                      "Important",
                                      "Critical"), 
                           ordered = TRUE)

d$Architects.involvement <- factor(d$Architects.involvement,
                                   levels = c("None",
                                              "Simple",
                                              "Monitoring",
                                              "Active Participation",
                                              "Joint Design"), 
                                   ordered = TRUE)

# binary
d$Critical.feature <- ifelse(d$Critical.feature == 'Yes', 1, 0)
d$Dependency <- ifelse(d$Dependency == 'Yes', 1, 0)

Descriptive statistics

We now have a data frame, d, which has all variable types correctly set.

str(d)
## 'data.frame':    11110 obs. of  10 variables:
##  $ ID                    : num  1 2 3 4 5 6 7 8 9 10 ...
##  $ State                 : Ord.factor w/ 6 levels "Elicited, Dropped"<..: 2 2 6 2 1 2 2 2 2 2 ...
##  $ Team.priority         : num  88 83 957 79 88 76 74 73 72 71 ...
##  $ Critical.feature      : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Business.value        : Ord.factor w/ 4 levels "No value"<"Valuable"<..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Customer.value        : Ord.factor w/ 4 levels "No value"<"Valuable"<..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Stakeholders          : num  1 1 0 0 0 1 1 1 1 1 ...
##  $ Key.customers         : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Dependency            : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Architects.involvement: Ord.factor w/ 5 levels "None"<"Simple"<..: 1 1 1 1 1 1 1 1 1 1 ...

Let’s plot our outcome and predictors so we get a feeling for the distributions. In the margin you will find all variables plotted.

First, we see that for State approximately as many features are released (final stage) as dropped in the first state. We also see that it drops off after the initial state.

For Team.priority many features have zero in priority (\(5139\)), and then there’s a bunch of them (\(1516\)) that have priority set to the maximum value, i.e., \(1000\).

For Critical.feature we have a clear emphasis on ‘No’.

Concerning Business.value and Customer.value they are fairly similar in their respective distribution (as one would expect).

For Stakeholder and Key.customers we see an emphasis on lower numbers, while for Dependency a clear emphasis on ‘No’.

Finally, for Architects.involvement we see that in the absolute majority of the cases architects are not involved.

In short, it looks sort of what one would expect, i.e., it’s not hard to find answers to why the plots look the way they do.

However, before we continue, we should standardize some of our predictors so the sampling will be easier, i.e., we simply do \((x - \bar{x})/\sigma_x\), then simply multiplying with \(\sigma_x\) and adding the mean, will allow us to get back to the original scale. It’s good practice to store this in new variables and suffix them with _s. At the same time, let’s give our variables shorter names.

# standardize and abbreviated names
d$prio_s <- scale(d$Team.priority)
d$sh_s <- scale(d$Stakeholders)
d$kc_s <- scale(d$Key.customers)

# abbreviate names
d$crit <- d$Critical.feature
d$b_val <- d$Business.value
d$c_val <- d$Customer.value
d$dep <- d$Dependency
d$arch <- d$Architects.involvement

Initial model comparison

First we sample a model that only has a population-level intercept (our null model, \(\mathcal{M}_0\)). Then we sample models using a cumulative likelihood (first w/o modeling predictors as monotonic) [2[2] P.-C. Bürkner and M. Vuorre, “Ordinal regression models in psychology: A tutorial,” Advances in Methods and Practices in Psychological Science, vol. 2, no. 1, pp. 77–101, 2019, doi: 10.1177/2515245918823199.]. Once we’ve sampled our models we then use LOO to compare each model’s relative out of sample prediction capabilities [3[3] G. Vehtari A. and J. Gabry, “Practical bayesian model evaluation using leave-one-out cross-validation and WAIC,” Statistics and Computing, vol. 27, pp. 1413–1432, 2017, doi: 10.1007/s11222-016-9696-4.]. Since our outcome is, in all practical sense, a Likert scale variable using an adjacent category (acat) model is really not an option, the acat model is a common ordinal model in item-response theory. Nor is a sequential model an option, i.e., we’re not after predicting the number of released features.

M0 <- brm(State ~ 1, family = cumulative, data = d, refresh = 0)
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Users/torkarr/Library/R/4.0/library/Rcpp/include/"  -I"/Users/torkarr/Library/R/4.0/library/RcppEigen/include/"  -I"/Users/torkarr/Library/R/4.0/library/RcppEigen/include/unsupported"  -I"/Users/torkarr/Library/R/4.0/library/BH/include" -I"/Users/torkarr/Library/R/4.0/library/StanHeaders/include/src/"  -I"/Users/torkarr/Library/R/4.0/library/StanHeaders/include/"  -I"/Users/torkarr/Library/R/4.0/library/rstan/include" -DEIGEN_NO_DEBUG  -D_REENTRANT  -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -include stan/math/prim/mat/fun/Eigen.hpp   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Users/torkarr/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/Core:88:
## /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Users/torkarr/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
## /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
M1 <- brm(State ~ 1 + prio_s + crit + b_val + c_val + sh_s + kc_s + dep + arch,
          family = cumulative, data = d, refresh = 0)
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Users/torkarr/Library/R/4.0/library/Rcpp/include/"  -I"/Users/torkarr/Library/R/4.0/library/RcppEigen/include/"  -I"/Users/torkarr/Library/R/4.0/library/RcppEigen/include/unsupported"  -I"/Users/torkarr/Library/R/4.0/library/BH/include" -I"/Users/torkarr/Library/R/4.0/library/StanHeaders/include/src/"  -I"/Users/torkarr/Library/R/4.0/library/StanHeaders/include/"  -I"/Users/torkarr/Library/R/4.0/library/rstan/include" -DEIGEN_NO_DEBUG  -D_REENTRANT  -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -include stan/math/prim/mat/fun/Eigen.hpp   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Users/torkarr/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/Core:88:
## /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Users/torkarr/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
## /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
# Since it's cumulative likelihood we could model as monotonic.
M2 <- brm(State ~ 1 + prio_s + crit + mo(b_val) + mo(c_val) + 
                 sh_s + kc_s + dep + mo(arch),
          family = cumulative, data = d, refresh = 0)
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Users/torkarr/Library/R/4.0/library/Rcpp/include/"  -I"/Users/torkarr/Library/R/4.0/library/RcppEigen/include/"  -I"/Users/torkarr/Library/R/4.0/library/RcppEigen/include/unsupported"  -I"/Users/torkarr/Library/R/4.0/library/BH/include" -I"/Users/torkarr/Library/R/4.0/library/StanHeaders/include/src/"  -I"/Users/torkarr/Library/R/4.0/library/StanHeaders/include/"  -I"/Users/torkarr/Library/R/4.0/library/rstan/include" -DEIGEN_NO_DEBUG  -D_REENTRANT  -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -include stan/math/prim/mat/fun/Eigen.hpp   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Users/torkarr/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/Core:88:
## /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Users/torkarr/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
## /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1

Compare the three models’ out of sample prediction capabilities.

(l <- loo_compare(loo(M0), loo(M1), loo(M2)))
##    elpd_diff se_diff
## M2     0.0       0.0
## M1    -2.1       2.9
## M0 -3170.2      62.1

LOO puts \(\mathcal{M}_2\) as no. 1. If we assume a \(z_{\text{95%}}\)-score of \(1.96\) it’s clear that zero is in the interval and that \(\mathcal{M}_2\) does not have much of an advantage, i.e., \(\text{CI}_{z_{95\%}}\) [-7.77, 3.57]. We can conclude this matter by saying that adding predictors to \(\mathcal{M}_0\) clearly has a significant effect, but adding monotonic effects to our model does not. Interpreting monotonic effects is slightly harder so in this case, since it has no significant effect, we will not include it in the model.

If we would be interested in refining our model for out of sample prediction purposes we could conduct variable selection. However, in this particular case we are interested in each predictor’s effect, so we’ll keep them all, and decide to use \(\mathcal{M}_1\) as our target model, \(\mathcal{M}\), for now.

Prior predictive checks

Before we start to use our data to do inferences we should think about our priors. Our outcome of interest is State, which is ordered categorical data and as such we will model it using a cumulative likelihood (as seen in the previous section).

Let’s see what priors we would get with a full model where we include all variables as predictors. We’ll model b_val, c_val, and arch as population-level effects, since the model comparison provided indications that there was no advantage to model these variables as monotonic effects.

We have 15 \(\beta\) and 5 \(\alpha\) parameters that we need to set a prior on. Since we use a \(\mathop{\mathrm{logit}}\) link for our model (to translate to the probabilistic scale) it is hard to conceptually grasp the effect of our priors, i.e., we need to visualize them.

First, we tried weakly regularizing priors on the \(\beta\) and \(\alpha\) parameters.1 Basically \(\mathcal{N}(0,2)\) and \(\mathcal{N}(0,5)\) Then, after analyzing and comparing different priors we conclude that the following priors are suitable for this model.2 Providing all steps in a prior sensitivity analysis is beyond the scope of this document and seldom reported in scientific literature.

p <- get_prior(State ~ 1 + prio_s + crit + b_val + c_val + sh_s + kc_s + dep + arch,
          family = cumulative, data = d)

# Beta (we have fifteen of them)
# M(0,0.1) might seem very strict, but (15*0.1)^2 = 2.25
p$prior[1] <- "normal(0, 0.1)"

# Our priors on alpha (five of them; one for each border between 
# Likert scale 1,...,6)
p$prior[17] <- "normal(0, 2)"

Let’s sample from the priors only (it’s enough to use one chain here), and then check against the data to see that they are wide enough (but not too wide!)

M_prior <- brm(State ~ 1 + prio_s + crit + b_val + c_val + sh_s + kc_s + dep + arch,
          family = cumulative, data = d, prior = p, sample_prior = "only", 
          chains = 1, refresh = 0)
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Users/torkarr/Library/R/4.0/library/Rcpp/include/"  -I"/Users/torkarr/Library/R/4.0/library/RcppEigen/include/"  -I"/Users/torkarr/Library/R/4.0/library/RcppEigen/include/unsupported"  -I"/Users/torkarr/Library/R/4.0/library/BH/include" -I"/Users/torkarr/Library/R/4.0/library/StanHeaders/include/src/"  -I"/Users/torkarr/Library/R/4.0/library/StanHeaders/include/"  -I"/Users/torkarr/Library/R/4.0/library/rstan/include" -DEIGEN_NO_DEBUG  -D_REENTRANT  -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -include stan/math/prim/mat/fun/Eigen.hpp   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Users/torkarr/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/Core:88:
## /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Users/torkarr/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
## /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1

pp_check(M_prior, type = "bars", nsamples = NULL)

If we sample from our prior predictive distribution, \(y_{\text{rep}}\), and plot it with our empirical data, \(y\), we’ll see how the priors cover our data. As is evident (see right), by looking at the intervals, our priors seems to be relaxed and should not influence the outcome, and in our case we have \(11100\) rows so priors should not affect our outcome either way. Additionally, the priors are nearly uniform on the outcome scale, which they should be; and the only sane way to check this is to plot the combination of all priors used in a model.

Inference

Let’s now sample from our model using the priors from above.

M <- brm(State ~ 1 + prio_s + crit + b_val + c_val + sh_s + kc_s + dep + arch, 
         family = cumulative, data = d, prior = p, chains = 4, refresh = 0)
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Users/torkarr/Library/R/4.0/library/Rcpp/include/"  -I"/Users/torkarr/Library/R/4.0/library/RcppEigen/include/"  -I"/Users/torkarr/Library/R/4.0/library/RcppEigen/include/unsupported"  -I"/Users/torkarr/Library/R/4.0/library/BH/include" -I"/Users/torkarr/Library/R/4.0/library/StanHeaders/include/src/"  -I"/Users/torkarr/Library/R/4.0/library/StanHeaders/include/"  -I"/Users/torkarr/Library/R/4.0/library/rstan/include" -DEIGEN_NO_DEBUG  -D_REENTRANT  -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -include stan/math/prim/mat/fun/Eigen.hpp   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Users/torkarr/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/Core:88:
## /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Users/torkarr/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
## /Users/torkarr/Library/R/4.0/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1

Posterior predictive check

Again, plotting our model’s predictions, \(y_{\text{rep}}\), vs. our empirical values, \(y\), we see that the model does a good job in estimating the categories. The credible intervals are tight (a lot of data). If we look at our previous plot we can clearly see, when we compare with the plot to the right, that our data has completely swamped the priors. We see some under- and overestimation in some categories but these are very minor.

For those so inclined, \(R^2 =\).

Diagnostics

Our caterpillar plots look good for all parameters the model estimated.

Our \(\widehat{R}\) and effective sample size (neff) look good.

# should be < 1.01
max(rhat(M), na.rm = TRUE)
## [1] 1.001304
# should be > 0.1
min(neff_ratio(M), na.rm = TRUE)
## [1] 0.4881814

There is no perfect model, but we could claim that this is a useful model. We would even consider taking it out for a dinner.

Estimates of parameters

If we now focus on our estimates we see that our ‘Intercepts’, which in this case are the borders between two categories in our outcome, are precisely estimated. However, to make sense of them we need to transform the estimates, since the values are on the \(\mathop{\mathrm{logit}}\) scale.

Our population-level estimates (fixed effects) from our model \(\mathcal{M}\). Note that the values are on \(\mathop{\mathrm{logit}}\) scale.

Estimate Est.Error Q2.5 Q97.5
Intercept[1] -1.68 0.08 -1.84 -1.53
Intercept[2] -0.46 0.07 -0.61 -0.32
Intercept[3] 0.76 0.07 0.61 0.90
Intercept[4] 1.15 0.07 1.00 1.29
Intercept[5] 1.20 0.07 1.06 1.35
prio_s 1.43 0.02 1.38 1.48
crit 0.66 0.05 0.56 0.76
b_val.L 0.20 0.05 0.11 0.30
c_val.L 0.11 0.06 0.00 0.22
sh_s -0.12 0.02 -0.16 -0.08
kc_s 0.02 0.02 -0.02 0.05
dep 0.11 0.05 0.02 0.20
arch.L 0.11 0.08 -0.05 0.25

Let’s take the first parameter, Intercept[1], which was estimated to -1.68, i.e.,

inv_logit_scaled(fixef(M)[1,1])
## [1] 0.1567248

What does 0.16 mean? Well, we have a ordered categorical outcome. So 16% of the probability mass, with a 95% credible interval of [0.14, 0.18], was assigned to the first category: Elicited, Dropped. For Intercept[2], 0.39 [0.35, 0.42] was assigned to: Elicited, Dropped and Elicited, Prio, Dropped.

Let’s turn our attention to the other parameters. The prio_s parameter, 1.43, would then become 0.81 when transformed. But remember the suffix _s! We need to multiply with \(\sigma_x\) and add \(\bar{x}\) from the data, which leads to an estimate of 704.29 [701.24, 707.24].

However, looking at point estimates is, quite frankly, not terribly useful. Let’s plot the posterior probability densitities for our population-level estimates on the \(\mathop{\mathrm{logit}}\) scale, disregarding Intercept\([1,\ldots,5]\).

Examining the above plot, from top to bottom, we can say that on the 95%-level, the first three parameters are clearly positive and do not cross \(0\). The fourth parameter, Customer value, is actually significant, with the following 95% credible intervals [0, 0.22]. The fifth parameter is clearly negative. The sixth parameter, Key customers, is not significant. Finally, Dependency is positive, while Architects involvement is not significant.

To conclude what we’ve noticed so far: Team priority, Critical feature, Business value, Stakeholders, and Dependency are clearly significant on \(\text{CI}_{95\%}\), while Customer value is borderline significant. Architecture involvement and Key customers are not significant.

Conditional effects

Below we plot all conditional effects for our model. The colors represent the different categories, \(1,\ldots,6\), for our outcome State. We are particularly interested in category 6 (pink), i.e., the final category which indicates a released feature.

  1. Elicited, Dropped
  2. Elicited, Prio, Dropped
  3. Elicited, Prio, Planned, Dropped
  4. Elicited, Prio, Planned, Implemented, Dropped
  5. Elicited, Prio, Planned, Implemented, Tested, Dropped
  6. Elicited, Prio, Planned, Implemented, Tested, Released

Take-aways concerning released features

One important question we would like to have an answer to is which independent variable(s) contribute(s) more for a feature to, ultimately, be released, i.e., is it priority, criticality, business or customer value, number of stakeholders, number of key customers, having dependencies, and/or the level of architect involvement? In the above plots the answer to our question can be found, without even having to conduct any statistical tests or examining \(p\)-values.

Concerning Priority we see that it has a very large effect for State 6 (i.e., a feature being released). The higher the priority the more probability mass is set on State 6. In the end it has close to 70% of the probability mass, while the other states are not even close.

Concerning Criticality we see very much the same effect, albeit the uncertainty increases also. States 3 and 6 have, together, more than 50% of the probability mass.

Business value and Customer value play a small role if we look at the plots above. However, one thing we clearly see in the Stakeholders plot is that State 6 has very little probability mass when number of stakeholders increase.

Key customers doesn’t tell us much due to the very large uncertainty in State 6. For Dependencies not much changes when it increases.

Finally, for Architects involvement we see that the probability for State 6 increases when going from none to Simple, but then it remains virtually the same for the other categories.

In short, Priority and Criticality have the largest effects, while the rest don’t matter all too much.